from scipy.linalg import solve
# Define the equation parameters
from scipy.special import erfc
import pandas as pd
import numpy as np
import scipy.integrate
from numpy import exp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.animation as FuncAnimation
from IPython.display import HTML
plt.rcParams['figure.figsize'] = [10, 5]
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
import seaborn as sns
from matplotlib.animation import PillowWriter

# Analytical - du/dt = -v*du/dx+D*d2u/dx2 - c*u +d

Ci = 0; 
Co = 1;
# dCinf/dx = 0;
R = 1; 
c = 0.0
D = 0.01
v = 1.0
d = 0.0
nx = 101;
nt = 1001;

L = 1.0;
T = 1.0;

x = np.linspace(0, L, nx);
t = np.linspace(0, T, nt);


C = np.zeros((nt, nx))

for i in range(nt):  
    for j in range(nx):
        C[i, j] = (0.025/(np.sqrt(0.000625 + 0.02*t[i])))*np.exp((-(x[j] + 0.5 - t[i])**2)/(0.00125 + 0.04*t[i]))
        
        # An = np.exp(-c*t[i]/R)*(1-(erfc((R*x[j]-v*t[i])/(2*np.sqrt(D*R*t[i]))))*0.5 -  0.5*(np.exp(v*x[j]/(D)))*(erfc((R*x[j]+v*t[i])/(2*np.sqrt(D*R*t[i])))))       
        # un = v*(np.sqrt(1+ (4*c*D)/(v**2)))
        # Bn = 0.5*np.exp(((v-un)*x[j])/(2*D))*(erfc((R*x[j]-un*t[i])/(2*np.sqrt(D*R*t[i])))) + 0.5*np.exp(((v+un)*x[j])/(2*D))*(erfc((R*x[j] + un*t[i])/(2*np.sqrt(D*R*t[i]))))
        # print(np.shape(An), np.shape(Bn))
        # C[i, j] = d/c + (Ci - d/c)*An + (Co-d/c)*Bn

U_exact = np.transpose(C);

fig = plt.figure()  
plt.imshow(np.transpose(C), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
plt.colorbar()
plt.title("Exact")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDI.png', dpi = 900)
plt.show()

[1.38389653e-87 1.45844014e-84 1.00924437e-81 ... 4.21713645e-04
 4.13846376e-04 4.06120826e-04]



[0.         0.         0.         ... 0.00038267 0.00039424 0.00040612]

from mpl_toolkits.mplot3d import Axes3D
sns.set(style = "darkgrid")
ax = fig.add_subplot(111, projection = '3d')
X, Y = np.meshgrid(x, np.transpose(t))
sns.set_style("whitegrid", {'axes.grid' : False})
fig = plt.figure(figsize=(12,12))
ax = plt.axes(projection='3d')
ax.scatter(X, Y, np.transpose(U_exact), s=4, c=U_exact, cmap = 'viridis', alpha = 1)
ax.set_xlabel('Position')
ax.set_ylabel('Time')
ax.set_zlabel('U_exact')
plt.show()

$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x}$$$$ 0<x<1 \ \ \ 0<t \leq 1 $$$$ IC: U(x, 0) = \exp(- \frac {(x+0.5)^2}{0.00125}) $$

$$ BC - Left: U(0, t) = \frac{0.025}{\sqrt{0.000625 + 0.02t}} \exp\left[ - \frac{(0.5-t)^2}{(0.00125 + 0.04t)} \right] $$ $$ BC - Right: U(1, t) = \frac{0.025}{\sqrt{0.000625 + 0.02t}} \exp \left[- \frac{(1.5-t)^2}{(0.00125 + 0.04t)} \right] $$ $$ 0<t \leq 1 $$

Analytical Solution - Para et al. 2022

For v = 1.0 and D = 0.01 $$ U(x, t) = \frac{0.025}{\sqrt{0.000625 + 0.02t}} \exp \left[- \frac{(x+0.5-t)^2}{(0.00125 + 0.04t)}\right] $$

# Define the grid
a = -1.0;
b = 0.01;
c = 0.0;
d = 0.0;



nx = 101
nt = 1001
L = 1.0
T = 1.0
# Define x & t
x = np.linspace(0, L, nx)
t = np.linspace(0, T, nt)


# Dirichlet Boundary Condition
left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)

fig = plt.figure()  
plt.imshow(np.transpose(C), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
plt.colorbar()
plt.title("Exact")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDI.png', dpi = 900)
plt.show()


# Discretization of space and time
dx = L/(nx-1)
dt = T/(nt-1)
u = np.exp(-((x+0.5)**2)/(0.00125)) #np.zeros(nx);
def numerical_scheme_CDI(u_CDI, dx, dt, l, r):
    alpha = a*dt/dx
    beta = b*dt/dx**2
    A = np.zeros((len(u_CDI), len(u_CDI)))
    bi = d*dt+u_CDI;

    for j in range(1, len(u_CDI)-1):

        A[j, j-1] = -(beta-alpha/2)
        A[j, j] = 1+2*beta-c*dt
        A[j, j+1] = -(beta+alpha/2)

    N = len(u_CDI)-1
    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0
    
    
    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r


    u_CDI = solve(A, bi)

    return u_CDI

def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r):

    frac = 0.1

    le = -np.log(1-l*frac)
    re = -np.log(1-r*frac)
    ke = -np.log(1-u_eCDI*frac)
    N = len(u_eCDI)-1
    #plt.plot(u_eCDI, label='u_eCDI')
    #plt.plot(ke, label='ke')
    A = np.zeros((len(u_eCDI), len(u_eCDI)))
    udterm = -c
    ydterm = d
    vdterm = -a
    wdterm = b
    #plt.legend()
    uk = udterm*dt*np.ones(N+1); # f
    u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
    vk = -vdterm*dt*np.ones(N+1); # a
    v2k = -wdterm*dt*np.ones(N+1); #d
    wk = wdterm*dt*np.ones(N+1); # b


    v2t = np.zeros(N+1)

    v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

    ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
    be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2));
    ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2)

    de = np.zeros(N+1);

    for t in range(0, (N+1)):
        if(t==N):
            de[t] = ke[t];
        else:
            de[t] = ke[t]*(-u2k[t]*np.exp(ke[t])/2+1)+uk[t]+u2k[t]*np.exp(ke[t])

    bi = de;

    for j in range(1, len(u_eCDI)-1):
        A[j, j-1] = ce[j]
        A[j, j] = 1+ae[j]
        A[j, j+1] = be[j]

    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0


    bi[0] = le
    bi[N] = re


    ke = solve(A, bi)
    u_eCDI = ((1-np.exp(-ke))/frac)
    return u_eCDI

def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r):

    # frac = 0.1


    # ke = -np.log(1-u_eCDCN*frac)
    # N = len(u_eCDCN)-1

    # A = np.zeros((len(u_eCDCN), len(u_eCDCN)))
    # udterm = -c
    # ydterm = d
    # vdterm = -a
    # wdterm = b


    # uk = udterm*dt*np.ones(N+1); # f
    # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
    # vk = -vdterm*dt*np.ones(N+1); # a
    # v2k = -wdterm*dt*np.ones(N+1); #d
    # wk = wdterm*dt*np.ones(N+1); # b


    # v2t = np.zeros(N+1)

    # v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

    # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
    # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2));
    # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2))

    # dcn = np.zeros(N+1);
    # for t in range(0, N+1):
    #     if(t==(N)):
    #         dcn[t] = ke[t];
    #     elif(t==0):
    #         dcn[t] = ke[t];
    #     else:
    #         dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


    frac = 0.1

    le = -np.log(1-l*frac)
    re = -np.log(1-r*frac)
    ke = -np.log(1-u_eCDCN*frac)
    N = len(u_eCDCN)-1

    A = np.zeros((len(u_eCDCN), len(u_eCDCN)))

    v2t = np.zeros(N+1)

    v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) 

    bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t)
    acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4)
    ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) -  ((b*dt)/(8*(dx**2)))*(v2t)

    dcn = np.zeros(N+1);
    for t in range(0, N+1):
        if(t==(N)):
            dcn[t] = ke[t];
        elif(t==0):
            dcn[t] = ke[t];
        else:
            dcn[t] = ke[t] - c*dt + ((d*frac + c)*dt*np.exp(ke[t]))*(1-ke[t]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[t+1] - ke[t-1]) + ((b*dt)/(dx**2))*(ke[t+1]+ke[t-1]-2*ke[t]) - ((b*dt)/(4*(dx**2)))*((ke[t+1] - ke[t-1])**2))
            #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


    bi = dcn;

    for j in range(1, N):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = ccn[j]
        A[j, j] = 1+acn[j]
        A[j, j+1] = bcn[j]

    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = le
    bi[N] = re

    # print(A)
    ke = solve(A, bi)
    u_eCDCN = ((1-np.exp(-ke))/frac)
    return u_eCDCN

def numerical_scheme_CDCN(u_CDCN, dx, dt, l, r):
    alpha = a*dt/(4*dx)
    beta = b*dt/(2*dx**2)

    A = np.zeros((len(u_CDCN), len(u_CDCN)))

    bi = np.zeros(len(u_CDCN))
    for j in range(1, len(u_CDCN)-1):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = -(beta - alpha)
        A[j, j] = 1+(2*beta - c*dt/2)
        A[j, j+1] = -(beta + alpha)

        bi[j] = d*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha)

    N = len(u)-1
    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r



    u_CDCN = solve(A, bi)

    return u_CDCN

def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r):
    frac = 0.1
    le = -np.log(1-l*frac)
    re = -np.log(1-r*frac)
    ke = -np.log(1-u_eQI*frac)
    N = len(u_eQI)-1
    #plt.plot(u_eCDI, label='u_eCDI')
    #plt.plot(ke, label='ke')
    A = np.zeros((N+1, N+1))
    udterm = -c
    ydterm = d
    vdterm = -a
    wdterm = b

    uk = udterm*dt*np.ones(N+1); # f
    u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
    vk = -vdterm*dt*np.ones(N+1); # a
    v2k = -wdterm*dt*np.ones(N+1); #d
    wk = wdterm*dt*np.ones(N+1); # b


    # yk = (bia*frac)*dt;#c+e+uk   *np.ones((1, int(N/2)+1)) #d
    v2t = np.zeros(N+1)
    # print(np.shape(v2t))
    v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2])

    v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3])
    # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)])
    # print(np.shape(v2t))
    #v2t = np.concatenate((v2t, v2t[0]))
    #print(np.shape(v2t))


    cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2)
    bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2)
    aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2)
    fq = vk/(8*dx) + v2k*v2t/((8*dx)**2)

    # USING LAST node Q with other side variation
    cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2)
    bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2)
    anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2)
    fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2)

    dq = np.zeros(N+1);
    for t in range(0, N+1):
        if(t==N):
            dq[t] = (ke[t])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2))
        else:
            dq[t] = ke[t]*((-u2k[t]*np.exp(ke[t])/2)+1)+uk[t]+u2k[t]*np.exp(ke[t]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2]


    bi = dq;

    for j in range(1, N-1):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = cq[j]
        A[j, j] = 1+aq[j]
        A[j, j+1] = bq[j]
        A[j, j+2] = fq[j]


    A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one
    A[N-1, N-1] = 1+anq[N-1]
    A[N-1, N] = bnq[N-1]
    A[N-1, N-3] = fnq[N-1]

    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = le
    bi[N] = re


    # print(A)
    ke = solve(A, bi)
    u_eQI = ((1-np.exp(-ke))/frac)
    return u_eQI

def numerical_scheme_QI(u_QI, dx, dt, l, r):
    alpha = a*dt
    beta = b*dt

    A = np.zeros((len(u_QI), len(u_QI)))
    bi = d*dt+u_QI;
    cQ = 3*alpha/(8*dx) - beta/(dx**2)
    bQ = -7*alpha/(8*dx) - beta/(dx**2)
    aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    eQ = 1*alpha/(8*dx)

    cnQ = 7*alpha/(8*dx) - beta/(dx**2)
    bnQ = -3*alpha/(8*dx) - beta/(dx**2)
    anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    enQ = -1*alpha/(8*dx)


    for j in range(1, len(u_QI)-2):

        A[j, j-1] = cQ
        A[j, j] = 1+aQ
        A[j, j+1] = bQ
        A[j, j+2] = eQ

    N = len(u_QI)-1

    A[N-1, N-2] = cnQ  # applied the a>0 terms for this one
    A[N-1, N-1] = 1+anQ
    A[N-1, N] = bnQ
    A[N-1, N-3] = enQ


    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r



    u_QI = solve(A, bi)

    return u_QI


def numerical_scheme_QCN(u_QCN, dx, dt, l, r):
    alpha = a*dt
    beta = b*dt

    A = np.zeros((len(u_QCN), len(u_QCN)))
    bi = d*dt+u_QCN;
    cQ = 3*alpha/(8*dx) - beta/(dx**2)
    bQ = -7*alpha/(8*dx) - beta/(dx**2)
    aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    eQ = 1*alpha/(8*dx)

    cnQ = 7*alpha/(8*dx) - beta/(dx**2)
    bnQ = -3*alpha/(8*dx) - beta/(dx**2)
    anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    enQ = -1*alpha/(8*dx)


    for j in range(1, len(u_QCN)-2):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = cQ/2
        A[j, j] = 1+aQ/2
        A[j, j+1] = bQ/2
        A[j, j+2] = eQ/2
        bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2

    N = len(u_QCN)-1

    A[N-1, N-2] = cnQ/2  # applied the a>0 terms for this one
    A[N-1, N-1] = 1+anQ/2
    A[N-1, N] = bnQ/2
    A[N-1, N-3] = enQ/2
    bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2


    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r



    u_QCN = solve(A, bi)

    return u_QCN


# u_init = u_exact(x, 1.0)

# Apply the numerical scheme

# U_FD = np.zeros((nx, nt))
U_CDI = np.zeros((nx, nt))
U_CDCN = np.zeros((nx, nt))
U_QI = np.zeros((nx, nt))
U_QCN = np.zeros((nx, nt))
U_eCDI = np.zeros((nx, nt))
U_eQI = np.zeros((nx, nt))

U_eCDCN = np.zeros((nx, nt))


# error_2norm_FD_CDI = np.zeros(nt)
error_2norm_CDI_CDCN = np.zeros(nt)
error_2norm_CDCN_QI = np.zeros(nt)
error_2norm_CDCN_eCDCN = np.zeros(nt)
error_2norm_CDCN_QI = np.zeros(nt)
error_2norm_QI_eQI = np.zeros(nt)

error_2norm_CDI_eCDI = np.zeros(nt)
error_2norm_CDI_QI = np.zeros(nt)
error_2norm_QI_QCN = np.zeros(nt)
error_2norm_CDCN_QCN = np.zeros(nt)
error_2norm_CDI_QCN = np.zeros(nt)


error_2norm_eCDCN_eCDI = np.zeros(nt)
error_2norm_eCDCN_eQI = np.zeros(nt)
error_2norm_eCDI_eQI = np.zeros(nt)

error_2norm_Qexact_QI = np.zeros(nt)
error_2norm_Qexact_eQI = np.zeros(nt)

error_2norm_Qexact_CDCN = np.zeros(nt)

error_2norm_Qexact_CDI = np.zeros(nt)
#error_2norm_FD_CDI = np.zeros(nt)


# u_FD = u;
u_CDI=u;
u_CDCN = u;
u_eQI = u;
u_eCDI = u;

u_QI=u;
u_QCN=u;
u_eCDCN=u;

for n in range(nt):
    # print(n)
    # u_FD = numerical_scheme(u_FD, dx, dt)
    # #print(u)
    # #print(str(n*dt)+" "+str(max(u)))
    # U_FD[:, n] = u_FD
    u_CDI = numerical_scheme_CDI(u_CDI, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_CDI[:, n] = u_CDI

    #error = u_FD - u_CDI
    #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)




    u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c, d, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_eCDI[:, n] = u_eCDI

    error = u_CDI - u_eCDI
    error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_CDCN = numerical_scheme_CDCN(u_CDCN, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_CDCN[:, n] = u_CDCN
    error = u_CDI - u_CDCN
    error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c, d, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_eCDCN[:, n] = u_eCDCN
    error = u_CDCN - u_eCDCN
    error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_QI = numerical_scheme_QI(u_QI, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_QI[:, n] = u_QI
    error = u_CDCN - u_QI
    error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_CDI - u_QI
    error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_eQI[:, n] = u_eQI

    error = u_QI - u_eQI
    error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_eCDI - u_eQI
    error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_eCDCN - u_eQI
    error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_eCDCN - u_eCDI
    error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_QCN = numerical_scheme_QCN(u_QCN, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_QCN[:, n] = u_QCN
    error = u_CDCN - u_QCN
    error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_CDI - u_QCN
    error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_QI - u_QCN
    error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    # Comparing with the exact solution
    error = U_exact[:, n] - u_QI
    error_2norm_Qexact_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


    error = U_exact[:, n] - u_CDCN
    error_2norm_Qexact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


    error = U_exact[:, n] - u_CDI
    error_2norm_Qexact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)
    
    error = U_exact[:, n] - u_eQI
    error_2norm_Qexact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    
plt.rcParams['font.size'] = 18
# fig = plt.figure()  
# plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
# plt.colorbar()
# plt.title("FD")
# plt.xlabel('Time')
# plt.ylabel('Position')

# fig.set_size_inches(30.,18.)
# plt.savefig('U_FD.png', dpi = 900)
# plt.show()

fig = plt.figure()  
plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("CDI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDI.png', dpi = 900)
plt.show()


fig = plt.figure()  
plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("eCDI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_eCDI.png', dpi = 900)
plt.show()


fig = plt.figure() 
plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("CDCN")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDCN.png', dpi = 900)
plt.show()


fig = plt.figure() 
plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("eCDCN")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDCN.png', dpi = 900)
plt.show()


fig = plt.figure() 
plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("eQI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_eQI.png', dpi = 900)
plt.show()

fig = plt.figure() 
plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("QI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_QI.png', dpi = 900)
plt.show()

fig = plt.figure() 
plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("QCN")
plt.xlabel('Time')
plt.ylabel('Position')

# fig.set_size_inches(30.,18.)
# plt.savefig('U_QCN.png', dpi = 900)
plt.show()



fig = plt.figure() 
#plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN')
#plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI')
#plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI')
#plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI')
#plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN')

#plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN')
#plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN')
#plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI')
#plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI')
plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI')
plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI')

plt.plot(error_2norm_Qexact_QI, '*', label = 'Exact_QI')


plt.plot(error_2norm_Qexact_CDI, '+', label = 'Exact_CDI')
plt.plot(error_2norm_Qexact_CDCN, 'x', label = 'Exact_CDDCN')

plt.plot(error_2norm_Qexact_eQI, '-', label = 'Exact_eQI')
plt.title("Errors")
plt.legend()
plt.xlabel('Iterations in Time')
plt.ylabel('Error-2-Norm')

#fig.set_size_inches(30.,18.)
# plt.savefig('Error_Norm.png', dpi = 900)
plt.show()

fig = plt.figure() 
# plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI')
#plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN')
#plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI')
#plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI')
#plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN')

#plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI')
#plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN')
#plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN')
#plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI')

#plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI')
#plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI')
#plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI')
plt.plot(np.log(error_2norm_Qexact_QI), '*', label = 'Exact_QI')
plt.plot(np.log(error_2norm_Qexact_CDI), '+', label = 'Exact_CDI')
plt.plot(np.log(error_2norm_Qexact_CDCN), 'x', label = 'Exact_CDCN')
plt.plot(np.log(error_2norm_Qexact_eQI), '-', label = 'Exact_eQI')

plt.legend()
plt.title("Log errors")
plt.xlabel('Iterations in Time')
plt.ylabel('log-Error-2-Norm')

#fig.set_size_inches(30.,18.)
# plt.savefig('Log-ErrorNorm.png', dpi = 900)
plt.show()


#fig = plt.figure() 
#plt.plot(U_exact[:, -1], '-+',label = 'U_exact')
#plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI')
#plt.legend()
#plt.title("Comparison of U_exact and U_CDI")
#plt.xlabel('Length')
#plt.ylabel('U')

#fig.set_size_inches(30.,18.)
# plt.savefig('Log-ErrorNorm.png', dpi = 900)
#plt.show()

plt.plot(U_eQI[:, -100], label='eQI'); plt.plot(U_exact[:, -100], label='Exact'); plt.legend(); plt.title("At 0.1 seconds - Pe = 100")

Text(0.5, 1.0, 'At 0.1 seconds - Pe = 100')

Adding MMS code with nx -> 20, 40, 80, 160, 240

For equation with IC and BCs:

$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x}$$$$ IC: U(x, 0) = C_i $$$$ BC - Left: U(0, t) = C_o $$$$ BC - Right: \frac {dU(1, t)}{dx} = 0 $$

Notes on MMS:

* Values should be within the domain of the code: TAU or here u/C lies in the range of 0-5 (0-1 for $$ Pe~1 $$) 
In this case, note that whatever manufactured solution is picked must be continuously differentiable through its second spatial derivative and first temporal derivative. In addition, for illustrative purposes, the manufactured solution is chosen so that time will appear in the boundary conditions. Something of this nature can be used as well: 
T(x, t) = Tg0+Tgxcos(ATgxπx)+Tgtcos(ATgtπt)+Tgtxcos(ATgtxπtx)
The chosen equation is: 

$$ U(x,t) = C_o \cos(\pi x) e^{-t}$$

For x and t from 0 to 1: A = 1 and B = pi/2

xx = np.linspace(0, 1, 101)
tt = np.linspace(0, 1, 1001)
uu = np.zeros((101, 1001))
Co = 2;
for i in range(1001):
    for j in range(101):
        uu[j, i] =  Co*np.cos(np.pi*xx[j])*np.exp(-tt[i])
        
        
plt.imshow(uu, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x210ac871c88>

Implementing that in the equation for the derivative terms are:

$$ \frac {dU}{dt} = -C_o\cos(\pi x) e^{-t}$$$$ \frac {dU}{dx} = -C_o\pi \sin(\pi x) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o\pi^2 \cos(\pi x) e^{-t}$$

The initial and boundary condition terms: Here left is Dirichlet and right is Neumann

$$ U(x, 0) = C_o\cos(\pi x)$$$$ U(0, t) = C_oe^{-t}$$$$ \frac {dU}{dx} (1, t) = -C_o\sin(\pi 1) e^{-t} = 0$$

The new source term will be:

$$ d = \frac {dU}{dt} + v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2}$$$$ d = -C_o\cos(\pi x) e^{-t} - v*C_o\pi \sin(\pi x) e^{-t} + D*C_o\pi^2 \cos(\pi x) e^{-t}$$$$ d = -C_o e^{-t}(\cos(\pi x) + v*\pi \sin(\pi x) - D*\pi^2 \cos(\pi x)) $$

Notes on MMS for transient flows As a first step, to examine temporal order, spatial discretization must be held constant (fixed grid); then, to examine spatial order, temporal discretization must be held constant (fixed time step). Holding either discretization at a constant level introduces a fixed error that applies to all calculations in that group.
The error in such cases is of the form:

$$ DE = DE_temporal + DE_spatial $$

For transient, keeping one of the errors constant, the other can be examined. However, the error we get is sum of constant term, k, with the error, let's say spatial. Hence, to get the the order, the observed order $ \hat{p} $ must be calculated to remove the constant term with three levels instead of the two required for steady flows as:

$$ \hat{p} = \frac {\ln( \frac{||DE_{l+2}|| - ||DE_{l+1}||}{||DE_{l+1}|| - ||DE_l||})}{\ln (r)} $$

nx		nt
25		1000
50		1000
100		1000
200		1000
400		1000
100		250
100		500
100		2000
100		4000		$$ d = -C_o e^{-t}(\cos(\pi x) + v*\pi \sin(\pi x) - D*\pi^2 \cos(\pi x)) $$ # Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. def numericalcodes(nx, nt): # Define the grid a = -1.0; b = 0.01; c = 0.0; Co = 1.0; L = 1.0 T = 1.0 dm = np.zeros((nt, nx)); U_exact = np.zeros((nt, nx)); # Define x & t x = np.linspace(0, L, nx) t = np.linspace(0, T, nt) for i in range(nt): for j in range (nx): dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - a*np.pi*np.sin(np.pi*x[j]) - b*(np.pi**2)*np.cos(np.pi*x[j])) # Dirichlet Boundary Condition left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) right = np.ones(nt)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) for i in range(nt): for j in range (nx): U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i]) fig = plt.figure() plt.imshow(np.transpose(U_exact), cmap='viridis', extent=[0, T, L, 0], aspect='auto') plt.colorbar() plt.title("Exact") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() # Discretization of space and time dx = L/(nx-1) dt = T/(nt-1) u = Co*np.cos(np.pi*x) #np.zeros(nx); def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt, l, r): alpha = a*dt/dx beta = b*dt/dx**2 A = np.zeros((len(u_CDI), len(u_CDI))) bi = d*dt+u_CDI; for j in range(1, len(u_CDI)-1): A[j, j-1] = -(beta-alpha/2) A[j, j] = 1+2*beta-c*dt A[j, j+1] = -(beta+alpha/2) N = len(u_CDI)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_CDI = solve(A, bi) return u_CDI def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDI*frac) N = len(u_eCDI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((len(u_eCDI), len(u_eCDI))) udterm = -c ydterm = d vdterm = -a wdterm = b #plt.legend() uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2)); ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2) de = np.zeros(N+1); for xt in range(0, (N+1)): if(xt==N): de[xt] = ke[xt]; else: de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) bi = de; for j in range(1, len(u_eCDI)-1): A[j, j-1] = ce[j] A[j, j] = 1+ae[j] A[j, j+1] = be[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re ke = solve(A, bi) u_eCDI = ((1-np.exp(-ke))/frac) return u_eCDI def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r): # frac = 0.1 # ke = -np.log(1-u_eCDCN*frac) # N = len(u_eCDCN)-1 # A = np.zeros((len(u_eCDCN), len(u_eCDCN))) # udterm = -c # ydterm = d # vdterm = -a # wdterm = b # uk = udterm*dt*np.ones(N+1); # f # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c # vk = -vdterm*dt*np.ones(N+1); # a # v2k = -wdterm*dt*np.ones(N+1); #d # wk = wdterm*dt*np.ones(N+1); # b # v2t = np.zeros(N+1) # v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2)); # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2)) # dcn = np.zeros(N+1); # for t in range(0, N+1): # if(t==(N)): # dcn[t] = ke[t]; # elif(t==0): # dcn[t] = ke[t]; # else: # dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDCN*frac) N = len(u_eCDCN)-1 A = np.zeros((len(u_eCDCN), len(u_eCDCN))) v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t) acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4) ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) - ((b*dt)/(8*(dx**2)))*(v2t) dcn = np.zeros(N+1); for xt in range(0, N+1): if(xt==(N)): dcn[xt] = ke[xt]; elif(xt==0): dcn[xt] = ke[xt]; else: dcn[xt] = ke[xt] - c*dt + ((d[xt]*frac + c)*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2)) #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) bi = dcn; for j in range(1, N): # A[j, j-2] = -d[j]*dt A[j, j-1] = ccn[j] A[j, j] = 1+acn[j] A[j, j+1] = bcn[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re # print(A) ke = solve(A, bi) u_eCDCN = ((1-np.exp(-ke))/frac) return u_eCDCN def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r): alpha = a*dt/(4*dx) beta = b*dt/(2*dx**2) A = np.zeros((len(u_CDCN), len(u_CDCN))) bi = np.zeros(len(u_CDCN)) for j in range(1, len(u_CDCN)-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = -(beta - alpha) A[j, j] = 1+(2*beta - c*dt/2) A[j, j+1] = -(beta + alpha) bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha) N = len(u)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_CDCN = solve(A, bi) return u_CDCN def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eQI*frac) N = len(u_eQI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((N+1, N+1)) udterm = -c ydterm = d vdterm = -a wdterm = b uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b # yk = (bia*frac)*dt;#c+e+uk *np.ones((1, int(N/2)+1)) #d v2t = np.zeros(N+1) # print(np.shape(v2t)) v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2]) v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3]) # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)]) # print(np.shape(v2t)) #v2t = np.concatenate((v2t, v2t[0])) #print(np.shape(v2t)) cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2) bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2) aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2) fq = vk/(8*dx) + v2k*v2t/((8*dx)**2) # USING LAST node Q with other side variation cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2) bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2) anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2) fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2) dq = np.zeros(N+1); for xt in range(0, N+1): if(xt==N): dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2)) else: dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2] bi = dq; for j in range(1, N-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = cq[j] A[j, j] = 1+aq[j] A[j, j+1] = bq[j] A[j, j+2] = fq[j] A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one A[N-1, N-1] = 1+anq[N-1] A[N-1, N] = bnq[N-1] A[N-1, N-3] = fnq[N-1] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re # print(A) ke = solve(A, bi) u_eQI = ((1-np.exp(-ke))/frac) return u_eQI def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QI), len(u_QI))) bi = d*dt+u_QI; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QI)-2): A[j, j-1] = cQ A[j, j] = 1+aQ A[j, j+1] = bQ A[j, j+2] = eQ N = len(u_QI)-1 A[N-1, N-2] = cnQ # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ A[N-1, N] = bnQ A[N-1, N-3] = enQ # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_QI = solve(A, bi) return u_QI def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QCN), len(u_QCN))) bi = d*dt+u_QCN; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QCN)-2): # A[j, j-2] = -d[j]*dt A[j, j-1] = cQ/2 A[j, j] = 1+aQ/2 A[j, j+1] = bQ/2 A[j, j+2] = eQ/2 bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2 N = len(u_QCN)-1 A[N-1, N-2] = cnQ/2 # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ/2 A[N-1, N] = bnQ/2 A[N-1, N-3] = enQ/2 bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_QCN = solve(A, bi) return u_QCN # u_init = u_exact(x, 1.0) # Apply the numerical scheme # U_FD = np.zeros((nx, nt)) U_CDI = np.zeros((nx, nt)) U_CDCN = np.zeros((nx, nt)) U_QI = np.zeros((nx, nt)) U_QCN = np.zeros((nx, nt)) U_eCDI = np.zeros((nx, nt)) U_eQI = np.zeros((nx, nt)) U_eCDCN = np.zeros((nx, nt)) # error_2norm_FD_CDI = np.zeros(nt) error_2norm_CDI_CDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_CDCN_eCDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_QI_eQI = np.zeros(nt) error_2norm_CDI_eCDI = np.zeros(nt) error_2norm_CDI_QI = np.zeros(nt) error_2norm_QI_QCN = np.zeros(nt) error_2norm_CDCN_QCN = np.zeros(nt) error_2norm_CDI_QCN = np.zeros(nt) error_2norm_eCDCN_eCDI = np.zeros(nt) error_2norm_eCDCN_eQI = np.zeros(nt) error_2norm_eCDI_eQI = np.zeros(nt) error_2norm_Exact_QI = np.zeros(nt) error_2norm_Exact_eQI = np.zeros(nt) error_2norm_Exact_eCDCN = np.zeros(nt) error_2norm_Exact_eCDI = np.zeros(nt) error_2norm_Exact_CDCN = np.zeros(nt) error_2norm_Exact_QCN = np.zeros(nt) error_2norm_Exact_CDI = np.zeros(nt) #error_2norm_FD_CDI = np.zeros(nt) # u_FD = u; u_CDI=u; u_CDCN = u; u_eQI = u; u_eCDI = u; u_QI=u; u_QCN=u; u_eCDCN=u; for n in range(nt): # print(n) # u_FD = numerical_scheme(u_FD, dx, dt) # #print(u) # #print(str(n*dt)+" "+str(max(u))) # U_FD[:, n] = u_FD u_CDI = numerical_scheme_CDI(u_CDI,a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDI[:, n] = u_CDI #error = u_FD - u_CDI #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDI[:, n] = u_eCDI error = u_CDI - u_eCDI error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_CDCN = numerical_scheme_CDCN(u_CDCN, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDCN[:, n] = u_CDCN error = u_CDI - u_CDCN error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDCN[:, n] = u_eCDCN error = u_CDCN - u_eCDCN error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_QI = numerical_scheme_QI(u_QI, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QI[:, n] = u_QI error = u_CDCN - u_QI error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QI error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eQI[:, n] = u_eQI error = u_QI - u_eQI error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDI - u_eQI error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eQI error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eCDI error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_QCN = numerical_scheme_QCN(u_QCN, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QCN[:, n] = u_QCN error = u_CDCN - u_QCN error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QCN error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_QI - u_QCN error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) # Comparing with the exact solution error = U_exact[n, :] - u_QI error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx) errorQI = error_2norm_Exact_QI[n] error = U_exact[n, :] - u_QCN error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorQCN = error_2norm_Exact_QCN[n] error = U_exact[n, :] - u_eCDCN error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDCN = error_2norm_Exact_eCDCN[n] error = U_exact[n, :] - u_eCDI error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDI = error_2norm_Exact_eCDI[n] error = U_exact[n, :] - u_CDCN error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDCN = error_2norm_Exact_CDCN[n] error = U_exact[n, :] - u_CDI error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDI = error_2norm_Exact_CDI[n] error = U_exact[n, :] - u_eQI error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreQI = error_2norm_Exact_eQI[n] # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150) errornorms = np.ones(7) errornorms[0] = errorCDI errornorms[1] = errorCDCN errornorms[2] = errorQI errornorms[3] = errorQCN errornorms[4] = erroreCDI errornorms[5] = erroreCDCN errornorms[6] = erroreQI plt.rcParams['font.size'] = 18 # fig = plt.figure() # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto') # plt.colorbar() # plt.title("FD") # plt.xlabel('Time') # plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_FD.png', dpi = 900) # plt.show() fig = plt.figure() plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eCDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eQI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eQI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_QI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QCN") plt.xlabel('Time') plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_QCN.png', dpi = 900) plt.show() fig = plt.figure() #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN') #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI') #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI') #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI') #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN') #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN') #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN') #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI') #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI') # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI') # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI') plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI') plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI') plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN') plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN') plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI') plt.title("Errors") plt.legend() plt.xlabel('Iterations in Time') plt.ylabel('Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Error_Norm.png', dpi = 900) plt.show() fig = plt.figure() # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI') #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN') #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI') #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI') #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN') #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI') #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN') #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN') #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI') #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI') #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI') #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI') plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI') plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI') plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN') plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN') plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI') plt.legend() plt.title("Log errors") plt.xlabel('Iterations in Time') plt.ylabel('log-Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) plt.show() #fig = plt.figure() #plt.plot(U_exact[:, -1], '-+',label = 'U_exact') #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI') #plt.legend() #plt.title("Comparison of U_exact and U_CDI") #plt.xlabel('Length') #plt.ylabel('U') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) #plt.show() return errornorms # Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. def numericalcodes_dirichlet(nx, nt): # Define the grid a = -1.0; b = 0.01; c = 0.0; Co = 2.0; L = 1.0 T = 1.0 dm = np.zeros((nt, nx)); U_exact = np.zeros((nt, nx)); # Define x & t x = np.linspace(0, L, nx) t = np.linspace(0, T, nt) for i in range(nt): for j in range (nx): dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - a*np.pi*np.sin(np.pi*x[j]) - b*(np.pi**2)*np.cos(np.pi*x[j])) # Dirichlet Boundary Condition left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) right = -Co * np.exp(-t)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) for i in range(nt): for j in range (nx): U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i]) fig = plt.figure() plt.imshow(np.transpose(U_exact), cmap='viridis', extent=[0, T, L, 0], aspect='auto') plt.colorbar() plt.title("Exact") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() # Discretization of space and time dx = L/(nx-1) dt = T/(nt-1) u = Co*np.cos(np.pi*x) #np.zeros(nx); def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt, l, r): alpha = a*dt/dx beta = b*dt/dx**2 A = np.zeros((len(u_CDI), len(u_CDI))) bi = d*dt+u_CDI; for j in range(1, len(u_CDI)-1): A[j, j-1] = -(beta-alpha/2) A[j, j] = 1+2*beta-c*dt A[j, j+1] = -(beta+alpha/2) N = len(u_CDI)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 # A[N, N] = 3 # A[N, N-1] = -4 # A[N, N-2] = 1 # bi[N] = 0 A[0,0] = 1.0 A[N,N] = 1.0 bi[0] = l bi[N] = r u_CDI = solve(A, bi) return u_CDI def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDI*frac) N = len(u_eCDI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((len(u_eCDI), len(u_eCDI))) udterm = -c ydterm = d vdterm = -a wdterm = b #plt.legend() uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2)); ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2) de = np.zeros(N+1); for xt in range(0, (N+1)): if(xt==N): de[xt] = ke[xt]; else: de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) bi = de; for j in range(1, len(u_eCDI)-1): A[j, j-1] = ce[j] A[j, j] = 1+ae[j] A[j, j+1] = be[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 # A[N, N] = 3 # A[N, N-1] = -4 # A[N, N-2] = 1 # bi[N] = 0 A[0,0] = 1.0 A[N,N] = 1.0 bi[0] = le bi[N] = re ke = solve(A, bi) u_eCDI = ((1-np.exp(-ke))/frac) return u_eCDI def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r): # frac = 0.1 # ke = -np.log(1-u_eCDCN*frac) # N = len(u_eCDCN)-1 # A = np.zeros((len(u_eCDCN), len(u_eCDCN))) # udterm = -c # ydterm = d # vdterm = -a # wdterm = b # uk = udterm*dt*np.ones(N+1); # f # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c # vk = -vdterm*dt*np.ones(N+1); # a # v2k = -wdterm*dt*np.ones(N+1); #d # wk = wdterm*dt*np.ones(N+1); # b # v2t = np.zeros(N+1) # v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2)); # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2)) # dcn = np.zeros(N+1); # for t in range(0, N+1): # if(t==(N)): # dcn[t] = ke[t]; # elif(t==0): # dcn[t] = ke[t]; # else: # dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDCN*frac) N = len(u_eCDCN)-1 A = np.zeros((len(u_eCDCN), len(u_eCDCN))) v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t) acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4) ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) - ((b*dt)/(8*(dx**2)))*(v2t) dcn = np.zeros(N+1); for xt in range(0, N+1): if(xt==(N)): dcn[xt] = ke[xt]; elif(xt==0): dcn[xt] = ke[xt]; else: dcn[xt] = ke[xt] - c*dt + ((d[xt]*frac + c)*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2)) #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) bi = dcn; for j in range(1, N): # A[j, j-2] = -d[j]*dt A[j, j-1] = ccn[j] A[j, j] = 1+acn[j] A[j, j+1] = bcn[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 # A[N, N] = 3 # A[N, N-1] = -4 # A[N, N-2] = 1 # bi[N] = 0 A[0,0] = 1.0 A[N,N] = 1.0 bi[0] = le bi[N] = re # print(A) ke = solve(A, bi) u_eCDCN = ((1-np.exp(-ke))/frac) return u_eCDCN def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r): alpha = a*dt/(4*dx) beta = b*dt/(2*dx**2) A = np.zeros((len(u_CDCN), len(u_CDCN))) bi = np.zeros(len(u_CDCN)) for j in range(1, len(u_CDCN)-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = -(beta - alpha) A[j, j] = 1+(2*beta - c*dt/2) A[j, j+1] = -(beta + alpha) bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha) N = len(u)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 # A[N, N] = 3 # A[N, N-1] = -4 # A[N, N-2] = 1 # bi[N] = 0 A[0,0] = 1.0 A[N,N] = 1.0 bi[0] = l bi[N] = r u_CDCN = solve(A, bi) return u_CDCN def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eQI*frac) N = len(u_eQI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((N+1, N+1)) udterm = -c ydterm = d vdterm = -a wdterm = b uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b # yk = (bia*frac)*dt;#c+e+uk *np.ones((1, int(N/2)+1)) #d v2t = np.zeros(N+1) # print(np.shape(v2t)) v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2]) v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3]) # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)]) # print(np.shape(v2t)) #v2t = np.concatenate((v2t, v2t[0])) #print(np.shape(v2t)) cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2) bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2) aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2) fq = vk/(8*dx) + v2k*v2t/((8*dx)**2) # USING LAST node Q with other side variation cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2) bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2) anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2) fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2) dq = np.zeros(N+1); for xt in range(0, N+1): if(xt==N): dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2)) else: dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2] bi = dq; for j in range(1, N-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = cq[j] A[j, j] = 1+aq[j] A[j, j+1] = bq[j] A[j, j+2] = fq[j] A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one A[N-1, N-1] = 1+anq[N-1] A[N-1, N] = bnq[N-1] A[N-1, N-3] = fnq[N-1] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 # A[N, N] = 3 # A[N, N-1] = -4 # A[N, N-2] = 1 # bi[N] = 0 A[0,0] = 1.0 A[N,N] = 1.0 bi[0] = le bi[N] = re # print(A) ke = solve(A, bi) u_eQI = ((1-np.exp(-ke))/frac) return u_eQI def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QI), len(u_QI))) bi = d*dt+u_QI; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QI)-2): A[j, j-1] = cQ A[j, j] = 1+aQ A[j, j+1] = bQ A[j, j+2] = eQ N = len(u_QI)-1 A[N-1, N-2] = cnQ # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ A[N-1, N] = bnQ A[N-1, N-3] = enQ # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 # A[N, N] = 3 # A[N, N-1] = -4 # A[N, N-2] = 1 # bi[N] = 0 A[0,0] = 1.0 A[N,N] = 1.0 bi[0] = l bi[N] = r u_QI = solve(A, bi) return u_QI def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QCN), len(u_QCN))) bi = d*dt+u_QCN; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QCN)-2): # A[j, j-2] = -d[j]*dt A[j, j-1] = cQ/2 A[j, j] = 1+aQ/2 A[j, j+1] = bQ/2 A[j, j+2] = eQ/2 bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2 N = len(u_QCN)-1 A[N-1, N-2] = cnQ/2 # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ/2 A[N-1, N] = bnQ/2 A[N-1, N-3] = enQ/2 bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 # A[N, N] = 3 # A[N, N-1] = -4 # A[N, N-2] = 1 # bi[N] = 0 A[0,0] = 1.0 A[N,N] = 1.0 bi[0] = l bi[N] = r u_QCN = solve(A, bi) return u_QCN # u_init = u_exact(x, 1.0) # Apply the numerical scheme # U_FD = np.zeros((nx, nt)) U_CDI = np.zeros((nx, nt)) U_CDCN = np.zeros((nx, nt)) U_QI = np.zeros((nx, nt)) U_QCN = np.zeros((nx, nt)) U_eCDI = np.zeros((nx, nt)) U_eQI = np.zeros((nx, nt)) U_eCDCN = np.zeros((nx, nt)) # error_2norm_FD_CDI = np.zeros(nt) error_2norm_CDI_CDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_CDCN_eCDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_QI_eQI = np.zeros(nt) error_2norm_CDI_eCDI = np.zeros(nt) error_2norm_CDI_QI = np.zeros(nt) error_2norm_QI_QCN = np.zeros(nt) error_2norm_CDCN_QCN = np.zeros(nt) error_2norm_CDI_QCN = np.zeros(nt) error_2norm_eCDCN_eCDI = np.zeros(nt) error_2norm_eCDCN_eQI = np.zeros(nt) error_2norm_eCDI_eQI = np.zeros(nt) error_2norm_Exact_QI = np.zeros(nt) error_2norm_Exact_eQI = np.zeros(nt) error_2norm_Exact_eCDCN = np.zeros(nt) error_2norm_Exact_eCDI = np.zeros(nt) error_2norm_Exact_CDCN = np.zeros(nt) error_2norm_Exact_QCN = np.zeros(nt) error_2norm_Exact_CDI = np.zeros(nt) #error_2norm_FD_CDI = np.zeros(nt) # u_FD = u; u_CDI=u; u_CDCN = u; u_eQI = u; u_eCDI = u; u_QI=u; u_QCN=u; u_eCDCN=u; for n in range(nt): # print(n) # u_FD = numerical_scheme(u_FD, dx, dt) # #print(u) # #print(str(n*dt)+" "+str(max(u))) # U_FD[:, n] = u_FD u_CDI = numerical_scheme_CDI(u_CDI,a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDI[:, n] = u_CDI #error = u_FD - u_CDI #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDI[:, n] = u_eCDI error = u_CDI - u_eCDI error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_CDCN = numerical_scheme_CDCN(u_CDCN, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDCN[:, n] = u_CDCN error = u_CDI - u_CDCN error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDCN[:, n] = u_eCDCN error = u_CDCN - u_eCDCN error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_QI = numerical_scheme_QI(u_QI, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QI[:, n] = u_QI error = u_CDCN - u_QI error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QI error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eQI[:, n] = u_eQI error = u_QI - u_eQI error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDI - u_eQI error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eQI error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eCDI error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_QCN = numerical_scheme_QCN(u_QCN, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QCN[:, n] = u_QCN error = u_CDCN - u_QCN error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QCN error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_QI - u_QCN error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) # Comparing with the exact solution error = U_exact[n, :] - u_QI error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx) errorQI = error_2norm_Exact_QI[n] error = U_exact[n, :] - u_QCN error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorQCN = error_2norm_Exact_QCN[n] error = U_exact[n, :] - u_eCDCN error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDCN = error_2norm_Exact_eCDCN[n] error = U_exact[n, :] - u_eCDI error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDI = error_2norm_Exact_eCDI[n] error = U_exact[n, :] - u_CDCN error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDCN = error_2norm_Exact_CDCN[n] error = U_exact[n, :] - u_CDI error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDI = error_2norm_Exact_CDI[n] error = U_exact[n, :] - u_eQI error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreQI = error_2norm_Exact_eQI[n] # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150) errornorms = np.ones(7) errornorms[0] = errorCDI errornorms[1] = errorCDCN errornorms[2] = errorQI errornorms[3] = errorQCN errornorms[4] = erroreCDI errornorms[5] = erroreCDCN errornorms[6] = erroreQI plt.rcParams['font.size'] = 18 # fig = plt.figure() # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto') # plt.colorbar() # plt.title("FD") # plt.xlabel('Time') # plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_FD.png', dpi = 900) # plt.show() fig = plt.figure() plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eCDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eQI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eQI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_QI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QCN") plt.xlabel('Time') plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_QCN.png', dpi = 900) plt.show() fig = plt.figure() #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN') #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI') #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI') #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI') #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN') #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN') #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN') #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI') #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI') # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI') # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI') plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI') plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI') plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN') plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN') plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI') plt.title("Errors") plt.legend() plt.xlabel('Iterations in Time') plt.ylabel('Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Error_Norm.png', dpi = 900) plt.show() fig = plt.figure() # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI') #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN') #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI') #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI') #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN') #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI') #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN') #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN') #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI') #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI') #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI') #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI') plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI') plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI') plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN') plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN') plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI') plt.legend() plt.title("Log errors") plt.xlabel('Iterations in Time') plt.ylabel('log-Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) plt.show() #fig = plt.figure() #plt.plot(U_exact[:, -1], '-+',label = 'U_exact') #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI') #plt.legend() #plt.title("Comparison of U_exact and U_CDI") #plt.xlabel('Length') #plt.ylabel('U') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) #plt.show() return errornorms # Spatial MMS tests errorsx = np.ones((9, 7)) errornorms = np.ones(7) nx = 25; nt = 100; nxs = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400] nts = [250, 500, 1000, 2000, 4000] for i in range(9): errornorms = numericalcodes(nx, nt) errorsx[i, 0] = errornorms[0] errorsx[i, 1] = errornorms[1] errorsx[i, 2] = errornorms[2] errorsx[i, 3] = errornorms[3] errorsx[i, 4] = errornorms[4] errorsx[i, 5] = errornorms[5] errorsx[i, 6] = errornorms[6] nx = nx*2; print(errorsx) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI') plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN') plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI') plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN') plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI') plt.plot(nxs, np.log(errorsx[:, 5]), '-v', label = 'eCDCN') plt.plot(nxs, np.log(errorsx[:, 6]), '--', label = 'eQI') plt.xlabel('Grid Numbers') plt.ylabel('Log Error (MSE)') plt.legend() plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) px = np.zeros((7, 7)) for g in range(7): for k in range(7): px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.plot(nxs[2:9], px[:, 0], '-', label = 'CDI') plt.plot(nxs[2:9], px[:, 1], '-.', label = 'CDCN') plt.plot(nxs[2:9], px[:, 2], '-+', label = 'QI') plt.plot(nxs[2:9], px[:, 3], '-*', label = 'QCN') plt.plot(nxs[2:9], px[:, 4], '-x', label = 'eCDI') plt.plot(nxs[2:9], px[:, 5], '-v', label = 'eCDCN') plt.plot(nxs[2:9], px[:, 6], '--', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Observed Order') plt.legend() [[0.00268041 0.00154615 0.00139582 0.00280235 0.00280723 0.00142768 0.00144015] [0.00137882 0.00274691 0.00116189 0.00307689 0.00142022 0.0026693 0.00118179] [0.00114741 0.00307596 0.00110675 0.00315649 0.00116595 0.00301502 0.00112009] [0.00109689 0.00315977 0.00108789 0.00317959 0.00110969 0.00310298 0.00109941] [0.00108228 0.0031819 0.00108013 0.00318681 0.00109353 0.00312614 0.00109107] [0.00107715 0.0031881 0.00107662 0.00318933 0.00108795 0.0031326 0.00108735] [0.00107508 0.00319 0.00107495 0.00319031 0.00108575 0.00313456 0.0010856 ] [0.00107417 0.00319065 0.00107414 0.00319073 0.00108478 0.00313522 0.00108475] [0.00107375 0.0031909 0.00107374 0.00319092 0.00108434 0.00313548 0.00108433]] <matplotlib.legend.Legend at 0x1895eda6910> # Temporal MMS tests errorst = np.ones((8, 7)) errornorms = np.ones(7) nx = 100; nt = 250; nxs = [25, 50, 100, 200, 400] nts = [50, 100, 200, 400, 800, 1600, 3200, 6400] for i in range(8): errornorms = numericalcodes(nx, nt) errorst[i, 0] = errornorms[0] errorst[i, 1] = errornorms[1] errorst[i, 2] = errornorms[2] errorst[i, 3] = errornorms[3] errorst[i, 4] = errornorms[4] errorst[i, 5] = errornorms[5] errorst[i, 6] = errornorms[6] nt = nt*2; print(errorst) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) # fig.set_size_inches(10.,7.) plt.plot(nts, np.log(errorst[:, 0]), '-', label = 'CDI') plt.plot(nts, np.log(errorst[:, 1]), '.-', label = 'CDCN') plt.plot(nts, np.log(errorst[:, 2]), '-+', label = 'QI') plt.plot(nts, np.log(errorst[:, 3]), '-*', label = 'QCN') plt.plot(nts, np.log(errorst[:, 4]), '-x', label = 'eCDI') plt.plot(nts, np.log(errorst[:, 5]), '--', label = 'eCDCN') plt.plot(nts, np.log(errorst[:, 6]), '-v', label = 'eQI') plt.legend(loc = 'best') plt.xlabel('Grid Numbers') plt.ylabel('Log Error') plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) pt = np.zeros((6, 7)) for g in range(6): for k in range(7): pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nts[g]/nts[g+1])) plt.plot(nts[1:7], pt[:, 0], '-', label = 'CDI') plt.plot(nts[1:7], pt[:, 1], '.-', label = 'CDCN') plt.plot(nts[1:7], pt[:, 2], '-+', label = 'QI') plt.plot(nts[1:7], pt[:, 3], '-*', label = 'QCN') plt.plot(nts[1:7], pt[:, 4], '-x', label = 'eCDI') plt.plot(nts[1:7], pt[:, 5], '-v', label = 'eCDCN') plt.plot(nts[1:7], pt[:, 6], '--', label = 'eQI') plt.xlabel('Grid Numbers') plt.ylabel('Observed Order') plt.legend() [[4.80797305e-04 1.15463686e-03 4.40796293e-04 1.23423173e-03 4.92591522e-04 1.12685979e-03 4.47041682e-04] [2.70027692e-04 5.23881067e-04 2.23849185e-04 6.02638800e-04 2.79732016e-04 5.07387787e-04 2.27857313e-04] [1.74819852e-04 2.11660470e-04 1.17467654e-04 2.88579416e-04 1.83282460e-04 2.01049458e-04 1.20428469e-04] [1.35616045e-04 6.30726926e-05 6.61680044e-05 1.32099918e-04 1.43018264e-04 5.72280516e-05 6.86376449e-05] [1.20546357e-04 4.42000304e-05 4.26769721e-05 5.43335537e-05 1.27101641e-04 4.92834621e-05 4.48618322e-05] [1.14690533e-04 7.43220864e-05 3.28202700e-05 1.71079475e-05 1.20701354e-04 8.01449362e-05 3.47624483e-05] [1.12262471e-04 9.19383565e-05 2.89506673e-05 1.04211480e-05 1.17967154e-04 9.75857228e-05 3.06933237e-05] [1.11183181e-04 1.01009111e-04 2.74193734e-05 1.73996971e-05 1.16726232e-04 1.06528305e-04 2.90304857e-05]] <ipython-input-11-e7ef7324d6f6>:40: RuntimeWarning: invalid value encountered in log pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nts[g]/nts[g+1])) <matplotlib.legend.Legend at 0x1895ee497c0> # Printing the spatial grid tests import pandas as pd zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6])) df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Spatial25.csv', index=False) # Printing the temporal grid tests zipped2 = list(zip(nts, errorst[:, 0], errorst[:, 1], errorst[:, 2], errorst[:, 3], errorst[:, 4], errorst[:, 5], errorst[:, 6], pt[:, 0], pt[:, 1], pt[:, 2], pt[:, 3], pt[:, 4], pt[:, 5], pt[:, 6])) df2 = pd.DataFrame(zipped2, columns=['Temporal Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Temporal250.csv', index=False) # Dirichlet both sides for spatial # Spatial MMS tests errorsx = np.ones((9, 7)) errornorms = np.ones(7) nx = 25; nt = 2000; nxs = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400] nts = [250, 500, 1000, 2000, 4000] for i in range(9): errornorms = numericalcodes_dirichlet(nx, nt) errorsx[i, 0] = errornorms[0] errorsx[i, 1] = errornorms[1] errorsx[i, 2] = errornorms[2] errorsx[i, 3] = errornorms[3] errorsx[i, 4] = errornorms[4] errorsx[i, 5] = errornorms[5] errorsx[i, 6] = errornorms[6] nx = nx*2; print(errorsx) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 12}) plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI') plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN') plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI') plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN') plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI') plt.plot(nxs, np.log(errorsx[:, 5]), '-v', label = 'eCDCN') plt.plot(nxs, np.log(errorsx[:, 6]), '--', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel(' Log Error (MSE)') plt.legend() plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 12}) px = np.zeros((7, 7)) for g in range(7): for k in range(7): px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI') plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN') plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI') plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN') plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI') plt.plot(nxs[1:8], px[:, 5], '-v', label = 'eCDCN') plt.plot(nxs[1:8], px[:, 6], '--', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Observed Order') plt.legend() # Printing the spatial grid tests import pandas as pd zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6])) df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Spatial25_Dirichlet.csv', index=False) --------------------------------------------------------------------------- MemoryError Traceback (most recent call last) <ipython-input-15-78b68e555be1> in <module> 9 nts = [250, 500, 1000, 2000, 4000] 10 for i in range(9): ---> 11 errornorms = numericalcodes_dirichlet(nx, nt) 12 errorsx[i, 0] = errornorms[0] 13 errorsx[i, 1] = errornorms[1] <ipython-input-4-d04c1cfcce3c> in numericalcodes_dirichlet(nx, nt) 586 587 --> 588 u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, dm[n, :], left[n], right[n]) 589 #print(u) 590 #print(str(n*dt)+" "+str(max(u))) <ipython-input-4-d04c1cfcce3c> in numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r) 365 366 # print(A) --> 367 ke = solve(A, bi) 368 u_eQI = ((1-np.exp(-ke))/frac) 369 return u_eQI ~\anaconda3\lib\site-packages\scipy\linalg\basic.py in solve(a, b, sym_pos, lower, overwrite_a, overwrite_b, debug, check_finite, assume_a, transposed) 205 norm = '1' 206 --> 207 anorm = lange(norm, a1) 208 209 # Generalized case 'gesv' MemoryError: Unable to allocate 1.22 MiB for an array with shape (400, 400) and data type float64 # Spatial MMS tests 2 errorsx = np.ones((9, 7)) errornorms = np.ones(7) nx = 20; nt = 1000; nxs = [20, 40, 80, 160, 320, 640, 1280, 2560, 5120] nts = [250, 500, 1000, 2000, 4000] for i in range(9): errornorms = numericalcodes(nx, nt) errorsx[i, 0] = errornorms[0] errorsx[i, 1] = errornorms[1] errorsx[i, 2] = errornorms[2] errorsx[i, 3] = errornorms[3] errorsx[i, 4] = errornorms[4] errorsx[i, 5] = errornorms[5] errorsx[i, 6] = errornorms[6] nx = nx*2; print(errorsx) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI') plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN') plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI') plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN') plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI') plt.plot(nxs, np.log(errorsx[:, 5]), '-v', label = 'eCDCN') plt.plot(nxs, np.log(errorsx[:, 6]), '--', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Log Error (MSE)') plt.legend() plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) px = np.zeros((7, 7)) for g in range(7): for k in range(7): px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI') plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN') plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI') plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN') plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI') plt.plot(nxs[1:8], px[:, 5], '-v', label = 'eCDCN') plt.plot(nxs[1:8], px[:, 6], '--', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Observed Order') plt.legend() # Printing the spatial grid tests zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6])) df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Spatial20.csv', index=False) The initial and boundary condition terms: I.C.: Cisin(pix) Boundary Conditions: Periodic BC $$ U(x, 0) = C_o\cos( \pi(2x -1))$$ $$ U(0, t) = -C_o $$ $$ U(1, t) = -C_o $$ $$ \frac {dU}{dx} (0, t) = -2\pi C_o\sin( \pi(2*0 -1)) e^{-t} = 0$$ $$ \frac {dU}{dx} (1, t) = -2\pi C_o\sin(\pi(2*1 -1)) e^{-t} = 0$$ The derivative terms are: $$ \frac {dU}{dt} = -C_o\cos(\pi (2*x-1)) e^{-t}$$$$ \frac {dU}{dx} = -C_o 2 \pi \sin(\pi (2*x-1)) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o 4 \pi^2 \cos(\pi (2*x-1)) e^{-t}$$ The new source term will be: $$ d = \frac {dU}{dt} + v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2}$$$$ d = -C_o\cos(\pi (2*x-1)) e^{-t} - v*2*C_o\pi \sin(\pi (2*x-1)) e^{-t} + D*C_o*4*\pi^2 \cos(\pi (2*x-1)) e^{-t}$$$$ d = -C_o e^{-t}(\cos(\pi (2*x-1)) + v*2*\pi \sin(\pi (2*x-1)) - D*4*\pi^2 \cos(\pi (2*x-1))) $$ For the equation with zero and first order decay terms: Other details remain same, only d will change as shown below: $$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x} + c*U + do $$$$ IC: U(x, 0) = C_i $$$$ BC - Left: U(0, t) = C_o $$$$ BC - Right: \frac {dU(1, t)}{dx} = 0 $$ $$ The \ manufactured \ solution \ for \ MMS \ are: $$ $$ U(x,t) = C_o \cos(\pi x) e^{-t}$$$$ \frac {dU}{dt} = -C_o\cos(\pi x) e^{-t}$$$$ \frac {dU}{dx} = -C_o\pi \sin(\pi x) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o\pi^2 \cos(\pi x) e^{-t}$$$$ The \ initial \ and \ boundary \ conditions \ for \ MMS \ are:$$$$ U(x, 0) = C_o\cos(\pi x)$$$$ U(0, t) = C_oe^{-t}$$$$ \frac {dU}{dx} (1, t) = -C_o\sin(\pi 1) e^{-t} = 0$$ $$ The \ new \ form \ of \ d \ is \ as: $$$$ d = \frac {dU}{dt} + v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2} - c*U - do $$$$ d = -C_o\cos(\pi x) e^{-t} - v*C_o\pi \sin(\pi x) e^{-t} + D*C_o\pi^2 \cos(\pi x) e^{-t} - c*(C_o \cos(\pi x) e^{-t}) - do $$$$ d = -C_o e^{-t}(\cos(\pi x) + v*\pi \sin(\pi x) - D*\pi^2 \cos(\pi x) + c*\cos(\pi x)) - do $$ # Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. def numericalcodes_complete(nx, nt): # Define the grid a = 0.1; b = 0.001; c = 0.1; Co = 2.0; L = 1.0 T = 1.0 do = 0.0 dm = np.zeros((nt, nx)); U_exact = np.zeros((nt, nx)); # Define x & t x = np.linspace(0, L, nx) t = np.linspace(0, T, nt) for i in range(nt): for j in range (nx): dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - a*np.pi*np.sin(np.pi*x[j]) - b*(np.pi**2)*np.cos(np.pi*x[j]) + c*np.cos(np.pi*x[j])) - do # Dirichlet Boundary Condition left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) right = np.ones(nt)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) for i in range(nt): for j in range (nx): U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i]) fig = plt.figure() plt.imshow(np.transpose(U_exact), cmap='viridis', extent=[0, T, L, 0], aspect='auto') plt.colorbar() plt.title("Exact") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() # Discretization of space and time dx = L/(nx-1) dt = T/(nt-1) u = Co*np.cos(np.pi*x) #np.zeros(nx); def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt, l, r): alpha = a*dt/dx beta = b*dt/dx**2 A = np.zeros((len(u_CDI), len(u_CDI))) bi = d*dt+u_CDI; for j in range(1, len(u_CDI)-1): A[j, j-1] = -(beta-alpha/2) - c*dt/2 A[j, j] = 1+2*beta A[j, j+1] = -(beta+alpha/2) - c*dt/2 N = len(u_CDI)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_CDI = solve(A, bi) return u_CDI def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDI*frac) N = len(u_eCDI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((len(u_eCDI), len(u_eCDI))) udterm = -c ydterm = d vdterm = -a wdterm = b #plt.legend() uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2)); ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2) de = np.zeros(N+1); for xt in range(0, (N+1)): if(xt==N): de[xt] = ke[xt]; else: de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) bi = de; for j in range(1, len(u_eCDI)-1): A[j, j-1] = ce[j] A[j, j] = 1+ae[j] A[j, j+1] = be[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re ke = solve(A, bi) u_eCDI = ((1-np.exp(-ke))/frac) return u_eCDI def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r): # frac = 0.1 # ke = -np.log(1-u_eCDCN*frac) # N = len(u_eCDCN)-1 # A = np.zeros((len(u_eCDCN), len(u_eCDCN))) # udterm = -c # ydterm = d # vdterm = -a # wdterm = b # uk = udterm*dt*np.ones(N+1); # f # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c # vk = -vdterm*dt*np.ones(N+1); # a # v2k = -wdterm*dt*np.ones(N+1); #d # wk = wdterm*dt*np.ones(N+1); # b # v2t = np.zeros(N+1) # v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2)); # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2)) # dcn = np.zeros(N+1); # for t in range(0, N+1): # if(t==(N)): # dcn[t] = ke[t]; # elif(t==0): # dcn[t] = ke[t]; # else: # dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDCN*frac) N = len(u_eCDCN)-1 A = np.zeros((len(u_eCDCN), len(u_eCDCN))) v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t) acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4) ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) - ((b*dt)/(8*(dx**2)))*(v2t) dcn = np.zeros(N+1); for xt in range(0, N+1): if(xt==(N)): dcn[xt] = ke[xt]; elif(xt==0): dcn[xt] = ke[xt]; else: dcn[xt] = ke[xt] - c*dt + ((d[xt]*frac + c)*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2)) #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) bi = dcn; for j in range(1, N): # A[j, j-2] = -d[j]*dt A[j, j-1] = ccn[j] A[j, j] = 1+acn[j] A[j, j+1] = bcn[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re # print(A) ke = solve(A, bi) u_eCDCN = ((1-np.exp(-ke))/frac) return u_eCDCN def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r): alpha = a*dt/(4*dx) beta = b*dt/(2*dx**2) A = np.zeros((len(u_CDCN), len(u_CDCN))) bi = np.zeros(len(u_CDCN)) for j in range(1, len(u_CDCN)-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = -(beta - alpha) A[j, j] = 1+(2*beta - c*dt/2) A[j, j+1] = -(beta + alpha) bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha) N = len(u)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_CDCN = solve(A, bi) return u_CDCN def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eQI*frac) N = len(u_eQI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((N+1, N+1)) udterm = -c ydterm = d vdterm = -a wdterm = b uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b # yk = (bia*frac)*dt;#c+e+uk *np.ones((1, int(N/2)+1)) #d v2t = np.zeros(N+1) # print(np.shape(v2t)) v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2]) v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3]) # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)]) # print(np.shape(v2t)) #v2t = np.concatenate((v2t, v2t[0])) #print(np.shape(v2t)) cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2) bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2) aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2) fq = vk/(8*dx) + v2k*v2t/((8*dx)**2) # USING LAST node Q with other side variation cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2) bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2) anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2) fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2) dq = np.zeros(N+1); for xt in range(0, N+1): if(xt==N): dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2)) else: dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2] bi = dq; for j in range(1, N-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = cq[j] A[j, j] = 1+aq[j] A[j, j+1] = bq[j] A[j, j+2] = fq[j] A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one A[N-1, N-1] = 1+anq[N-1] A[N-1, N] = bnq[N-1] A[N-1, N-3] = fnq[N-1] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re # print(A) ke = solve(A, bi) u_eQI = ((1-np.exp(-ke))/frac) return u_eQI def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QI), len(u_QI))) bi = d*dt+u_QI; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QI)-2): A[j, j-1] = cQ A[j, j] = 1+aQ A[j, j+1] = bQ A[j, j+2] = eQ N = len(u_QI)-1 A[N-1, N-2] = cnQ # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ A[N-1, N] = bnQ A[N-1, N-3] = enQ # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_QI = solve(A, bi) return u_QI def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QCN), len(u_QCN))) bi = d*dt+u_QCN; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QCN)-2): # A[j, j-2] = -d[j]*dt A[j, j-1] = cQ/2 A[j, j] = 1+aQ/2 A[j, j+1] = bQ/2 A[j, j+2] = eQ/2 bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2 N = len(u_QCN)-1 A[N-1, N-2] = cnQ/2 # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ/2 A[N-1, N] = bnQ/2 A[N-1, N-3] = enQ/2 bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_QCN = solve(A, bi) return u_QCN # u_init = u_exact(x, 1.0) # Apply the numerical scheme # U_FD = np.zeros((nx, nt)) U_CDI = np.zeros((nx, nt)) U_CDCN = np.zeros((nx, nt)) U_QI = np.zeros((nx, nt)) U_QCN = np.zeros((nx, nt)) U_eCDI = np.zeros((nx, nt)) U_eQI = np.zeros((nx, nt)) U_eCDCN = np.zeros((nx, nt)) # error_2norm_FD_CDI = np.zeros(nt) error_2norm_CDI_CDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_CDCN_eCDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_QI_eQI = np.zeros(nt) error_2norm_CDI_eCDI = np.zeros(nt) error_2norm_CDI_QI = np.zeros(nt) error_2norm_QI_QCN = np.zeros(nt) error_2norm_CDCN_QCN = np.zeros(nt) error_2norm_CDI_QCN = np.zeros(nt) error_2norm_eCDCN_eCDI = np.zeros(nt) error_2norm_eCDCN_eQI = np.zeros(nt) error_2norm_eCDI_eQI = np.zeros(nt) error_2norm_Exact_QI = np.zeros(nt) error_2norm_Exact_eQI = np.zeros(nt) error_2norm_Exact_eCDCN = np.zeros(nt) error_2norm_Exact_eCDI = np.zeros(nt) error_2norm_Exact_CDCN = np.zeros(nt) error_2norm_Exact_QCN = np.zeros(nt) error_2norm_Exact_CDI = np.zeros(nt) #error_2norm_FD_CDI = np.zeros(nt) # u_FD = u; u_CDI=u; u_CDCN = u; u_eQI = u; u_eCDI = u; u_QI=u; u_QCN=u; u_eCDCN=u; for n in range(nt): # print(n) # u_FD = numerical_scheme(u_FD, dx, dt) # #print(u) # #print(str(n*dt)+" "+str(max(u))) # U_FD[:, n] = u_FD u_CDI = numerical_scheme_CDI(u_CDI,a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDI[:, n] = u_CDI #error = u_FD - u_CDI #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDI[:, n] = u_eCDI error = u_CDI - u_eCDI error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_CDCN = numerical_scheme_CDCN(u_CDCN, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDCN[:, n] = u_CDCN error = u_CDI - u_CDCN error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDCN[:, n] = u_eCDCN error = u_CDCN - u_eCDCN error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_QI = numerical_scheme_QI(u_QI, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QI[:, n] = u_QI error = u_CDCN - u_QI error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QI error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eQI[:, n] = u_eQI error = u_QI - u_eQI error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDI - u_eQI error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eQI error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eCDI error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_QCN = numerical_scheme_QCN(u_QCN, a, b, c, dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QCN[:, n] = u_QCN error = u_CDCN - u_QCN error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QCN error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_QI - u_QCN error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) # Comparing with the exact solution error = U_exact[n, :] - u_QI error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx) errorQI = error_2norm_Exact_QI[n] error = U_exact[n, :] - u_QCN error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorQCN = error_2norm_Exact_QCN[n] error = U_exact[n, :] - u_eCDCN error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDCN = error_2norm_Exact_eCDCN[n] error = U_exact[n, :] - u_eCDI error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDI = error_2norm_Exact_eCDI[n] error = U_exact[n, :] - u_CDCN error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDCN = error_2norm_Exact_CDCN[n] error = U_exact[n, :] - u_CDI error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDI = error_2norm_Exact_CDI[n] error = U_exact[n, :] - u_eQI error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreQI = error_2norm_Exact_eQI[n] # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150) errornorms = np.ones(7) errornorms[0] = errorCDI errornorms[1] = errorCDCN errornorms[2] = errorQI errornorms[3] = errorQCN errornorms[4] = erroreCDI errornorms[5] = erroreCDCN errornorms[6] = erroreQI plt.rcParams['font.size'] = 18 # fig = plt.figure() # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto') # plt.colorbar() # plt.title("FD") # plt.xlabel('Time') # plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_FD.png', dpi = 900) # plt.show() fig = plt.figure() plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eCDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eQI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eQI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_QI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QCN") plt.xlabel('Time') plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_QCN.png', dpi = 900) plt.show() fig = plt.figure() #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN') #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI') #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI') #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI') #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN') #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN') #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN') #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI') #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI') # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI') # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI') plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI') plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI') plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN') plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN') plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI') plt.title("Errors") plt.legend() plt.xlabel('Iterations in Time') plt.ylabel('Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Error_Norm.png', dpi = 900) plt.show() fig = plt.figure() # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI') #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN') #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI') #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI') #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN') #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI') #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN') #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN') #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI') #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI') #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI') #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI') plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI') plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI') plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN') plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN') plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI') plt.legend() plt.title("Log errors") plt.xlabel('Iterations in Time') plt.ylabel('log-Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) plt.show() #fig = plt.figure() #plt.plot(U_exact[:, -1], '-+',label = 'U_exact') #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI') #plt.legend() #plt.title("Comparison of U_exact and U_CDI") #plt.xlabel('Length') #plt.ylabel('U') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) #plt.show() plt.figure() # print(np.shape(U_exact)) # print(np.shape(U_eQI)) plt.plot(U_eQI[:, nt-1], label='eQI'); plt.plot(U_exact[nt-1, :], label='Exact'); plt.legend(); plt.title("At 0.1 seconds - Pe = 100") return errornorms # Equation with decay terms # Spatial MMS tests errorsx = np.ones((9, 7)) errornorms = np.ones(7) nx = 25; nt = 1000; nxs = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400] nts = [250, 500, 1000, 2000, 4000] for i in range(9): errornorms = numericalcodes_complete(nx+1, nt+1) errorsx[i, 0] = errornorms[0] errorsx[i, 1] = errornorms[1] errorsx[i, 2] = errornorms[2] errorsx[i, 3] = errornorms[3] errorsx[i, 4] = errornorms[4] errorsx[i, 5] = errornorms[5] errorsx[i, 6] = errornorms[6] nx = nx*2; print(errors) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) plt.plot(nxs, errorsx[:, 0], '-', label = 'CDI') plt.plot(nxs, errorsx[:, 1], '.-', label = 'CDCN') plt.plot(nxs, errorsx[:, 2], '-+', label = 'QI') plt.plot(nxs, errorsx[:, 3], '-*', label = 'QCN') plt.plot(nxs, errorsx[:, 4], '-x', label = 'eCDI') plt.plot(nxs, errorsx[:, 5], '-v', label = 'eCDCN') plt.plot(nxs, errorsx[:, 6], '--', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Log Error (MSE)') plt.legend() plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) px = np.zeros((7, 7)) for g in range(7): for k in range(7): px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.plot(nxs[1:6], px[:, 0], '-', label = 'CDI') plt.plot(nxs[1:6], px[:, 1], '-.', label = 'CDCN') plt.plot(nxs[1:6], px[:, 2], '-+', label = 'QI') plt.plot(nxs[1:6], px[:, 3], '-*', label = 'QCN') plt.plot(nxs[1:6], px[:, 4], '-x', label = 'eCDI') plt.plot(nxs[1:6], px[:, 5], '-v', label = 'eCDCN') plt.plot(nxs[1:6], px[:, 6], '--', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Observed Order') plt.legend() # Printing the spatial grid tests zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6])) df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Spatial25_Complete.csv', index=False) # Temporal MMS tests errorst = np.ones((8, 7)) errornorms = np.ones(7) nx = 100; nt = 50; nxs = [25, 50, 100, 200, 400] nts = [50, 100, 200, 400, 800, 1600, 3200, 6400] for i in range(5): errornorms = numericalcodes_complete(nx, nt) errorst[i, 0] = errornorms[0] errorst[i, 1] = errornorms[1] errorst[i, 2] = errornorms[2] errorst[i, 3] = errornorms[3] errorst[i, 4] = errornorms[4] errorst[i, 5] = errornorms[5] errorst[i, 6] = errornorms[6] nt = nt*2; print(errors) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) # fig.set_size_inches(10.,7.) plt.plot(nts, np.log(errorst[:, 0]), '-', label = 'CDI') plt.plot(nts, np.log(errorst[:, 1]), '.-', label = 'CDCN') plt.plot(nts, np.log(errorst[:, 2]), '-+', label = 'QI') plt.plot(nts, np.log(errorst[:, 3]), '-*', label = 'QCN') plt.plot(nts, np.log(errorst[:, 4]), '-x', label = 'eCDI') plt.plot(nts, np.log(errorst[:, 5]), '--', label = 'eCDCN') plt.plot(nts, np.log(errorst[:, 6]), '-v', label = 'eQI') plt.legend(loc = 'best') plt.xlabel('Number of Time-steps') plt.ylabel('Log Error (MSE)') plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) pt = np.zeros((6, 7)) for g in range(6): for k in range(7): pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nts[g]/nts[g+1])) plt.plot(nts[1:5], pt[:, 0], '-', label = 'CDI') plt.plot(nts[1:5], pt[:, 1], '.-', label = 'CDCN') plt.plot(nts[1:5], pt[:, 2], '-+', label = 'QI') plt.plot(nts[1:5], pt[:, 3], '-*', label = 'QCN') plt.plot(nts[1:5], pt[:, 4], '-x', label = 'eCDI') plt.plot(nts[1:5], pt[:, 5], '-v', label = 'eCDCN') plt.plot(nts[1:5], pt[:, 6], '--', label = 'eQI') plt.xlabel('Grid Numbers') plt.ylabel('Observed Order') plt.legend() # Printing the temporal grid tests zipped2 = list(zip(nts, errorst[:, 0], errorst[:, 1], errorst[:, 2], errorst[:, 3], errorst[:, 4], errorst[:, 5], errorst[:, 6], pt[:, 0], pt[:, 1], pt[:, 2], pt[:, 3], pt[:, 4], pt[:, 5], pt[:, 6])) df2 = pd.DataFrame(zipped2, columns=['Temporal Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Temporal250_complete.csv', index=False) # For ordered grid increases # Equation with decay terms # Spatial MMS tests errorsx = np.ones((9, 7)) errornorms = np.ones(7) nx = 20; nt = 1000; nxs = [20, 40, 80, 160, 320, 640, 1280, 2560, 5120] # nxs = [25, 50, 100, 200, 400] nts = [250, 500, 1000, 2000, 4000] for i in range(10): nx = nxs[i]; errornorms = numericalcodes_complete(int(nx)+1, nt+1) errorsx[i, 0] = errornorms[0] errorsx[i, 1] = errornorms[1] errorsx[i, 2] = errornorms[2] errorsx[i, 3] = errornorms[3] errorsx[i, 4] = errornorms[4] errorsx[i, 5] = errornorms[5] errorsx[i, 6] = errornorms[6] print(errors) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 12}) # fig.set_size_inches(10.,7.) plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI') plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN') plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI') plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN') plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI') plt.plot(nxs, np.log(errorsx[:, 5]), '--', label = 'eCDCN') plt.plot(nxs, np.log(errorsx[:, 6]), '-v', label = 'eQI') plt.legend(loc = 'best') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Log Error (MSE)') plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 12}) px = np.zeros((7, 7)) for g in range(7): for k in range(7): px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI') plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN') plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI') plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN') plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI') plt.plot(nxs[1:8], px[:, 5], '--', label = 'eCDCN') plt.plot(nxs[1:8], px[:, 6], '-v', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Observed Order') plt.legend(loc = 'best') # Printing the spatial grid tests zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6])) df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Spatial25_Complete_Ordered.csv', index=False) px = np.zeros((8, 7)) for g in range(8): for k in range(7): # print((np.log(nxs[g]/nxs[g+1]))) px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 12}) plt.plot(nxs[0:8], px[:, 0], '-', label = 'CDI') plt.plot(nxs[0:8], px[:, 1], '-.', label = 'CDCN') plt.plot(nxs[0:8], px[:, 2], '-+', label = 'QI') plt.plot(nxs[0:8], px[:, 3], '-*', label = 'QCN') plt.plot(nxs[0:8], px[:, 4], '-x', label = 'eCDI') plt.plot(nxs[0:8], px[:, 6], '-v', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Observed Order') plt.legend() np.log(0.5) Variable Coefficients a, b, c $$ a = \sin(\pi x)$$$$ b = \cos(\pi x) $$$$ c = (a+b)/2$$ # Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. def numericalcodes_variable(nx, nt): # Define the grid Co = 2.0; L = 1.0 T = 1.0 do = 0.0 dm = np.zeros((nt, nx)); U_exact = np.zeros((nt, nx)); # Define x & t x = np.linspace(0, L, nx) t = np.linspace(0, T, nt) am = np.zeros((nt, nx)) bm = np.zeros((nt, nx)) cm = np.zeros((nt, nx)) for i in range(nt): for j in range (nx): am[i, j] = 0.1+0.1*np.sin(np.pi*x[j])# *np.exp(-t[i]); bm[i, j] = 0.001+0.001*np.cos(np.pi*x[j])# *np.exp(-t[i]); cm[i, j] = (am[i, j]+bm[i, j])/2; dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - am[i, j]*np.pi*np.sin(np.pi*x[j]) - bm[i, j]*(np.pi**2)*np.cos(np.pi*x[j]) + cm[i, j]*np.cos(np.pi*x[j])) - do # Dirichlet Boundary Condition left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) right = np.ones(nt)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt) for i in range(nt): for j in range (nx): U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i]) fig = plt.figure() plt.imshow(np.transpose(U_exact), cmap='viridis', extent=[0, T, L, 0], aspect='auto') plt.colorbar() plt.title("Exact") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() # Discretization of space and time dx = L/(nx-1) dt = T/(nt-1) u = Co*np.cos(np.pi*x) #np.zeros(nx); def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt, l, r): # print(a,b,c,d,l) alpha = a*dt/dx beta = b*dt/dx**2 A = np.zeros((len(u_CDI), len(u_CDI))) bi = d*dt+u_CDI; for j in range(1, len(u_CDI)-1): A[j, j-1] = -(beta[j]-alpha[j]/2) - c[j]*dt/2 A[j, j] = 1+2*beta[j] A[j, j+1] = -(beta[j]+alpha[j]/2) - c[j]*dt/2 N = len(u_CDI)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_CDI = solve(A, bi) return u_CDI def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDI*frac) N = len(u_eCDI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((len(u_eCDI), len(u_eCDI))) udterm = -c ydterm = d vdterm = -a wdterm = b #plt.legend() uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2)); ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2) de = np.zeros(N+1); for xt in range(0, (N+1)): if(xt==N): de[xt] = ke[xt]; else: de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) bi = de; for j in range(1, len(u_eCDI)-1): A[j, j-1] = ce[j] A[j, j] = 1+ae[j] A[j, j+1] = be[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re ke = solve(A, bi) u_eCDI = ((1-np.exp(-ke))/frac) return u_eCDI def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r): # frac = 0.1 # ke = -np.log(1-u_eCDCN*frac) # N = len(u_eCDCN)-1 # A = np.zeros((len(u_eCDCN), len(u_eCDCN))) # udterm = -c # ydterm = d # vdterm = -a # wdterm = b # uk = udterm*dt*np.ones(N+1); # f # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c # vk = -vdterm*dt*np.ones(N+1); # a # v2k = -wdterm*dt*np.ones(N+1); #d # wk = wdterm*dt*np.ones(N+1); # b # v2t = np.zeros(N+1) # v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2)); # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2)) # dcn = np.zeros(N+1); # for t in range(0, N+1): # if(t==(N)): # dcn[t] = ke[t]; # elif(t==0): # dcn[t] = ke[t]; # else: # dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eCDCN*frac) N = len(u_eCDCN)-1 A = np.zeros((len(u_eCDCN), len(u_eCDCN))) v2t = np.zeros(N+1) v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t) acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4) ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) - ((b*dt)/(8*(dx**2)))*(v2t) dcn = np.zeros(N+1); for xt in range(0, N+1): if(xt==(N)): dcn[xt] = ke[xt]; elif(xt==0): dcn[xt] = ke[xt]; else: dcn[xt] = ke[xt] - c[xt]*dt + ((d[xt]*frac + c[xt])*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a[xt]*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b[xt]*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b[xt]*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2)) #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2)) bi = dcn; for j in range(1, N): # A[j, j-2] = -d[j]*dt A[j, j-1] = ccn[j] A[j, j] = 1+acn[j] A[j, j+1] = bcn[j] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re # print(A) ke = solve(A, bi) u_eCDCN = ((1-np.exp(-ke))/frac) return u_eCDCN def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r): alpha = a*dt/(4*dx) beta = b*dt/(2*dx**2) A = np.zeros((len(u_CDCN), len(u_CDCN))) bi = np.zeros(len(u_CDCN)) for j in range(1, len(u_CDCN)-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = -(beta[j] - alpha[j]) A[j, j] = 1+(2*beta[j] - c[j]*dt/2) A[j, j+1] = -(beta[j] + alpha[j]) bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta[j] - c[j]*dt/2))+u_CDCN[j+1]*(alpha[j]+beta[j])+u_CDCN[j-1]*(beta[j]-alpha[j]) N = len(u)-1 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_CDCN = solve(A, bi) return u_CDCN def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r): frac = 0.1 le = -np.log(1-l*frac) re = -np.log(1-r*frac) ke = -np.log(1-u_eQI*frac) N = len(u_eQI)-1 #plt.plot(u_eCDI, label='u_eCDI') #plt.plot(ke, label='ke') A = np.zeros((N+1, N+1)) udterm = -c ydterm = d vdterm = -a wdterm = b uk = udterm*dt*np.ones(N+1); # f u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c vk = -vdterm*dt*np.ones(N+1); # a v2k = -wdterm*dt*np.ones(N+1); #d wk = wdterm*dt*np.ones(N+1); # b # yk = (bia*frac)*dt;#c+e+uk *np.ones((1, int(N/2)+1)) #d v2t = np.zeros(N+1) # print(np.shape(v2t)) v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2]) v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3]) # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)]) # print(np.shape(v2t)) #v2t = np.concatenate((v2t, v2t[0])) #print(np.shape(v2t)) cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2) bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2) aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2) fq = vk/(8*dx) + v2k*v2t/((8*dx)**2) # USING LAST node Q with other side variation cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2) bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2) anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2) fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2) dq = np.zeros(N+1); for xt in range(0, N+1): if(xt==N): dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2)) else: dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2] bi = dq; for j in range(1, N-1): # A[j, j-2] = -d[j]*dt A[j, j-1] = cq[j] A[j, j] = 1+aq[j] A[j, j+1] = bq[j] A[j, j+2] = fq[j] A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one A[N-1, N-1] = 1+anq[N-1] A[N-1, N] = bnq[N-1] A[N-1, N-3] = fnq[N-1] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = le # bi[N] = re # print(A) ke = solve(A, bi) u_eQI = ((1-np.exp(-ke))/frac) return u_eQI def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QI), len(u_QI))) bi = d*dt+u_QI; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QI)-2): A[j, j-1] = cQ[j] A[j, j] = 1+aQ[j] A[j, j+1] = bQ[j] A[j, j+2] = eQ[j] N = len(u_QI)-1 A[N-1, N-2] = cnQ[N-1] # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ[N-1] A[N-1, N] = bnQ[N-1] A[N-1, N-3] = enQ[N-1] # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_QI = solve(A, bi) return u_QI def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r): alpha = a*dt beta = b*dt A = np.zeros((len(u_QCN), len(u_QCN))) bi = d*dt+u_QCN; cQ = 3*alpha/(8*dx) - beta/(dx**2) bQ = -7*alpha/(8*dx) - beta/(dx**2) aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt eQ = 1*alpha/(8*dx) cnQ = 7*alpha/(8*dx) - beta/(dx**2) bnQ = -3*alpha/(8*dx) - beta/(dx**2) anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt enQ = -1*alpha/(8*dx) for j in range(1, len(u_QCN)-2): # A[j, j-2] = -d[j]*dt A[j, j-1] = cQ[j]/2 A[j, j] = 1+aQ[j]/2 A[j, j+1] = bQ[j]/2 A[j, j+2] = eQ[j]/2 bi[j] = bi[j] - cQ[j]*u_QCN[j-1]/2 - aQ[j]*u_QCN[j]/2 - bQ[j]*u_QCN[j+1]/2 - eQ[j]*u_QCN[j+2]/2 N = len(u_QCN)-1 A[N-1, N-2] = cnQ[N-1]/2 # applied the a>0 terms for this one A[N-1, N-1] = 1+anQ[N-1]/2 A[N-1, N] = bnQ[N-1]/2 A[N-1, N-3] = enQ[N-1]/2 bi[N-1] = bi[N-1] - cnQ[N-1]*u_QCN[N-2]/2 - anQ[N-1]*u_QCN[N-1]/2 - bnQ[N-1]*u_QCN[N]/2 - enQ[N-1]*u_QCN[N-3]/2 # Boundary condition # A[0, 0] = -3 # A[0, 1] = 4 # A[0, 2] = -1 # bi[0] = 0 A[N, N] = 3 A[N, N-1] = -4 A[N, N-2] = 1 bi[N] = 0 A[0,0] = 1.0 # A[N,N] = 1.0 bi[0] = l # bi[N] = r u_QCN = solve(A, bi) return u_QCN # u_init = u_exact(x, 1.0) # Apply the numerical scheme # U_FD = np.zeros((nx, nt)) U_CDI = np.zeros((nx, nt)) U_CDCN = np.zeros((nx, nt)) U_QI = np.zeros((nx, nt)) U_QCN = np.zeros((nx, nt)) U_eCDI = np.zeros((nx, nt)) U_eQI = np.zeros((nx, nt)) U_eCDCN = np.zeros((nx, nt)) # error_2norm_FD_CDI = np.zeros(nt) error_2norm_CDI_CDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_CDCN_eCDCN = np.zeros(nt) error_2norm_CDCN_QI = np.zeros(nt) error_2norm_QI_eQI = np.zeros(nt) error_2norm_CDI_eCDI = np.zeros(nt) error_2norm_CDI_QI = np.zeros(nt) error_2norm_QI_QCN = np.zeros(nt) error_2norm_CDCN_QCN = np.zeros(nt) error_2norm_CDI_QCN = np.zeros(nt) error_2norm_eCDCN_eCDI = np.zeros(nt) error_2norm_eCDCN_eQI = np.zeros(nt) error_2norm_eCDI_eQI = np.zeros(nt) error_2norm_Exact_QI = np.zeros(nt) error_2norm_Exact_eQI = np.zeros(nt) error_2norm_Exact_eCDCN = np.zeros(nt) error_2norm_Exact_eCDI = np.zeros(nt) error_2norm_Exact_CDCN = np.zeros(nt) error_2norm_Exact_QCN = np.zeros(nt) error_2norm_Exact_CDI = np.zeros(nt) #error_2norm_FD_CDI = np.zeros(nt) # u_FD = u; u_CDI=u; u_CDCN = u; u_eQI = u; u_eCDI = u; u_QI=u; u_QCN=u; u_eCDCN=u; for n in range(nt): # print(n) # u_FD = numerical_scheme(u_FD, dx, dt) # #print(u) # #print(str(n*dt)+" "+str(max(u))) # U_FD[:, n] = u_FD u_CDI = numerical_scheme_CDI(u_CDI,am[n, :], bm[n, :], cm[n, :], dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDI[:, n] = u_CDI #error = u_FD - u_CDI #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, am[n, :], bm[n, :], cm[n, :], dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDI[:, n] = u_eCDI error = u_CDI - u_eCDI error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_CDCN = numerical_scheme_CDCN(u_CDCN, am[n, :], bm[n, :], cm[n, :], dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_CDCN[:, n] = u_CDCN error = u_CDI - u_CDCN error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, am[n, :], bm[n, :], cm[n, :], dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eCDCN[:, n] = u_eCDCN error = u_CDCN - u_eCDCN error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) u_QI = numerical_scheme_QI(u_QI, am[n, :], bm[n, :], cm[n, :], dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QI[:, n] = u_QI error = u_CDCN - u_QI error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QI error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx) u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, am[n, :], bm[n, :], cm[n, :], dm[n, :], left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_eQI[:, n] = u_eQI error = u_QI - u_eQI error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDI - u_eQI error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eQI error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_eCDCN - u_eCDI error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) u_QCN = numerical_scheme_QCN(u_QCN, am[n, :], bm[n, :], cm[n, :], dm[n, :], dx, dt, left[n], right[n]) #print(u) #print(str(n*dt)+" "+str(max(u))) U_QCN[:, n] = u_QCN error = u_CDCN - u_QCN error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_CDI - u_QCN error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) error = u_QI - u_QCN error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) # Comparing with the exact solution error = U_exact[n, :] - u_QI error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx) errorQI = error_2norm_Exact_QI[n] error = U_exact[n, :] - u_QCN error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorQCN = error_2norm_Exact_QCN[n] error = U_exact[n, :] - u_eCDCN error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDCN = error_2norm_Exact_eCDCN[n] error = U_exact[n, :] - u_eCDI error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreCDI = error_2norm_Exact_eCDI[n] error = U_exact[n, :] - u_CDCN error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDCN = error_2norm_Exact_CDCN[n] error = U_exact[n, :] - u_CDI error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx) errorCDI = error_2norm_Exact_CDI[n] error = U_exact[n, :] - u_eQI error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx) erroreQI = error_2norm_Exact_eQI[n] # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150) errornorms = np.ones(7) errornorms[0] = errorCDI errornorms[1] = errorCDCN errornorms[2] = errorQI errornorms[3] = errorQCN errornorms[4] = erroreCDI errornorms[5] = erroreCDCN errornorms[6] = erroreQI plt.rcParams['font.size'] = 18 # fig = plt.figure() # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto') # plt.colorbar() # plt.title("FD") # plt.xlabel('Time') # plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_FD.png', dpi = 900) # plt.show() fig = plt.figure() plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eCDI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("CDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eCDCN") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_CDCN.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("eQI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_eQI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QI") plt.xlabel('Time') plt.ylabel('Position') #fig.set_size_inches(30.,18.) # plt.savefig('U_QI.png', dpi = 900) plt.show() fig = plt.figure() plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto') plt.colorbar() plt.title("QCN") plt.xlabel('Time') plt.ylabel('Position') # fig.set_size_inches(30.,18.) # plt.savefig('U_QCN.png', dpi = 900) plt.show() fig = plt.figure() #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN') #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI') #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI') #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI') #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN') #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN') #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN') #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI') #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI') # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI') # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI') plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI') plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI') plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN') plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN') plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI') plt.title("Errors") plt.legend() plt.xlabel('Iterations in Time') plt.ylabel('Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Error_Norm.png', dpi = 900) plt.show() fig = plt.figure() # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI') #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN') #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI') #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI') #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN') #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI') #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN') #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN') #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI') #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI') #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI') #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI') plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI') plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI') plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN') plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN') plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI') plt.legend() plt.title("Log errors") plt.xlabel('Iterations in Time') plt.ylabel('log-Error-2-Norm') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) plt.show() #fig = plt.figure() #plt.plot(U_exact[:, -1], '-+',label = 'U_exact') #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI') #plt.legend() #plt.title("Comparison of U_exact and U_CDI") #plt.xlabel('Length') #plt.ylabel('U') #fig.set_size_inches(30.,18.) # plt.savefig('Log-ErrorNorm.png', dpi = 900) #plt.show() plt.figure() # print(np.shape(U_exact)) # print(np.shape(U_eQI)) plt.plot(U_eQI[:, nt-1], label='eQI'); plt.plot(U_exact[nt-1, :], label='Exact'); plt.legend(); plt.title("At 0.1 seconds - Pe = 100") return errornorms # For ordered grid increases with variable coefficients with x and t # Equation with first order decay term # Spatial MMS tests errorsx = np.ones((9, 7)) errornorms = np.ones(7) nx = 20; nt = 1000; #nxs = np.linspace(20, 400, 10) nxs = [20, 40, 80, 160, 320, 640, 1280, 2560, 5120] nts = [250, 500, 1000, 2000, 4000] for i in range(10): nx = nxs[i] errornorms = numericalcodes_variable(int(nx)+1, nt+1) errorsx[i, 0] = errornorms[0] errorsx[i, 1] = errornorms[1] errorsx[i, 2] = errornorms[2] errorsx[i, 3] = errornorms[3] errorsx[i, 4] = errornorms[4] errorsx[i, 5] = errornorms[5] errorsx[i, 6] = errornorms[6] plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 12}) # fig.set_size_inches(10.,7.) plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI') plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN') plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI') plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN') plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI') plt.plot(nxs, np.log(errorsx[:, 5]), '--', label = 'eCDCN') plt.plot(nxs, np.log(errorsx[:, 6]), '-v', label = 'eQI') plt.legend(loc = 'best') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Log Error (MSE)') plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 12}) px = np.zeros((7, 7)) for g in range(7): for k in range(7): px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI') plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN') plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI') plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN') plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI') plt.plot(nxs[1:8], px[:, 5], '--', label = 'eCDCN') plt.plot(nxs[1:8], px[:, 6], '-v', label = 'eQI') plt.xlabel('Spatial Grid Numbers') plt.ylabel('Observed Order') plt.legend(loc = 'best') # Printing the spatial grid tests zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6])) df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Spatial25_Variable_Ordered.csv', index=False) # Temporal MMS tests errorst = np.ones((9, 7)) errornorms = np.ones(7) nx = 100; nt = 50; nxs = [25, 50, 100, 200, 400] nts = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400] for i in range(7): nt = nts[i] errornorms = numericalcodes_complete(nx, nt) errorst[i, 0] = errornorms[0] errorst[i, 1] = errornorms[1] errorst[i, 2] = errornorms[2] errorst[i, 3] = errornorms[3] errorst[i, 4] = errornorms[4] errorst[i, 5] = errornorms[5] errorst[i, 6] = errornorms[6] print(errors) plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) # fig.set_size_inches(10.,7.) plt.plot(nts, np.log(errorst[:, 0]), '-', label = 'CDI') plt.plot(nts, np.log(errorst[:, 1]), '.-', label = 'CDCN') plt.plot(nts, np.log(errorst[:, 2]), '-+', label = 'QI') plt.plot(nts, np.log(errorst[:, 3]), '-*', label = 'QCN') plt.plot(nts, np.log(errorst[:, 4]), '-x', label = 'eCDI') plt.plot(nts, np.log(errorst[:, 5]), '--', label = 'eCDCN') plt.plot(nts, np.log(errorst[:, 6]), '-v', label = 'eQI') plt.legend(loc = 'best') plt.xlabel('Grid Numbers') plt.ylabel('Log Error(MSE)') plt.figure() plt.rcParams["font.family"] = "Times New Roman" plt.rcParams.update({'font.size': 15}) # fig.set_size_inches(10.,7.) pt = np.zeros((7, 7)) for g in range(7): for k in range(7): pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nxs[g]/nxs[g+1])) plt.plot(nts[1:8], pt[:, 0], '-', label = 'CDI') plt.plot(nts[1:8], pt[:, 1], '.-', label = 'CDCN') plt.plot(nts[1:8], pt[:, 2], '-+', label = 'QI') plt.plot(nts[1:8], pt[:, 3], '-*', label = 'QCN') plt.plot(nts[1:8], pt[:, 4], '-x', label = 'eCDI') plt.plot(nts[1:8], pt[:, 5], '-v', label = 'eCDCN') plt.plot(nts[1:8], pt[:, 6], '--', label = 'eQI') plt.xlabel('Number of Time-Steps') plt.ylabel('Observed Order') plt.legend(loc = 'best') # Printing the temporal grid tests zipped2 = list(zip(nts, errorst[:, 0], errorst[:, 1], errorst[:, 2], errorst[:, 3], errorst[:, 4], errorst[:, 5], errorst[:, 6], pt[:, 0], pt[:, 1], pt[:, 2], pt[:, 3], pt[:, 4], pt[:, 5], pt[:, 6])) df2 = pd.DataFrame(zipped2, columns=['Temporal Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order']) df2.to_csv('Analysis_Temporal250_complete.csv', index=False) [[4.80797305e-04 1.15463686e-03 4.40796293e-04 1.23423173e-03 4.92591522e-04 1.12685979e-03 4.47041682e-04] [2.70027692e-04 5.23881067e-04 2.23849185e-04 6.02638800e-04 2.79732016e-04 5.07387787e-04 2.27857313e-04] [1.74819852e-04 2.11660470e-04 1.17467654e-04 2.88579416e-04 1.83282460e-04 2.01049458e-04 1.20428469e-04] [1.35616045e-04 6.30726926e-05 6.61680044e-05 1.32099918e-04 1.43018264e-04 5.72280516e-05 6.86376449e-05] [1.20546357e-04 4.42000304e-05 4.26769721e-05 5.43335537e-05 1.27101641e-04 4.92834621e-05 4.48618322e-05]] --------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-409-9df30d52c836> in <module> 34 for g in range(5): 35 for k in range(7): ---> 36 pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nxs[g]/nxs[g+1])) 37 38 plt.plot(nxs[0:5], pt[:, 0], '-', label = 'CDI') IndexError: list index out of range <Figure size 432x288 with 0 Axes>